{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 什么是” 最小二乘法” ?\n",
    "\n",
    "      定义：最小二乘法（又称最小平方法）是一种数学优化技术，它通过最小化误差的平方和寻找数据的最佳函数匹配。\n",
    "\n",
    "      作用：利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。\n",
    "\n",
    "      原则：以” 残差平方和最小” 确定直线位置 (在数理统计中，残差是指实际观察值与估计值之间的差)\n",
    "\n",
    "      数学公式：\n",
    "\n",
    " ![image1](1.png)\n",
    "\n",
    "      基本思路：对于一元线性回归模型, 假设从总体中获取了 n 组观察值（X1，Y1），（X2，Y2）， …，（Xn，Yn），对于平面中的这 n 个点，可以使用无数条曲线来拟合。而线性回归就是要求样本回归函数尽可能好地拟合这组值，也就是说，这条直线应该尽可能的处于样本数据的中心位置。因此，选择最佳拟合曲线的标准可以确定为：使总的拟合误差（即总残差）达到最小。 求最小，那就通过 对参数分别求导数联立方程组来解。\n",
    "\n",
    "最小二乘法是直接利用最小化误差平法和，来对参数求导，求得参数解，属于比较确定的值\n",
    "而 梯度下降法 ，属于迭代法，知道梯度 下降的方向刘，挨个去迭代。\n",
    "### 1. 最小二乘法的原理与要解决的问题　\n",
    "![image2](2.png)\n",
    "\n",
    "### 2. 最小二乘法的代数法解法\n",
    "最小二乘法是线性回归问题最常用的解决方法。我们这里使用到的是普通最小二乘法，英文为：Ordinary Least Squares，简称 OLS。最小二乘法中的「二乘」代表平方，最小二乘也就是最小平方。而这里的平方就是指代平方损失函数。首先，假设我们需要求解的是一元线性回归问题，及其对应的线性方程如下：\n",
    "\n",
    "![image3](equation--10-.svg)\n",
    "\n",
    "简单来讲，最小二乘法也就是求解平方损失函数最小值的方法。使用到高等数学中的知识推导如下。\n",
    "\n",
    "我们所知，平方损失函数的公式为：\n",
    "\n",
    "![image3](equation--7-.svg)\n",
    "\n",
    "由于最终的目标是求取平方损失函数 min(f)min(f) 最小时对应的 ww。所以，先求 ff 的一阶偏导数：\n",
    "\n",
    "![image3](equation--8-.svg)\n",
    "\n",
    "然后，我们令$ \\frac{\\partial f}{\\partial w_{0}}=0 $ 以及$ \\frac{\\partial f}{\\partial w_{1}}=0 $\n",
    "解得：\n",
    "\n",
    "![image3](equation--9-.svg)\n",
    "\n",
    "至此，我们就用代数的方式求解出线性方程的参数了。其中，$ w_0 $经常习惯称之为截距项。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们依据上面的推导过程，使用 Python 进行实现。这里，假设进行线性拟合的数据如下：\n",
    "\n",
    "x = [55, 71, 68, 87, 101, 87, 75, 78, 93, 73]\n",
    "y = [91, 101, 87, 109, 129, 98, 95, 101, 104, 93]\n",
    "根据公式实现参数计算函数 ols_algebra(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = [55, 71, 68, 87, 101, 87, 75, 78, 93, 73]\n",
    "y = [91, 101, 87, 109, 129, 98, 95, 101, 104, 93]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44.25604341391219, 0.7175629008386778)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#最小二乘法代数求解实现\n",
    "\n",
    "def ols_algebra(x, y):\n",
    "    n = len(x)\n",
    "    w1 = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x*x) - sum(x)*sum(x))\n",
    "    w0 = (sum(x*x)*sum(y) - sum(x)*sum(x*y))/(n*sum(x*x)-sum(x)*sum(x))\n",
    "    return w0, w1\n",
    "\n",
    "ols_algebra(np.array(x), np.array(y))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最小二乘法矩阵求解推导\n",
    "除了上面的代数推导过程。通常在实际情况中，更多会使用矩阵的方式来求线性拟合参数。因为在大规模运算中，矩阵的计算效率更高。\n",
    "\n",
    "矩阵的推导会稍微复杂一些。同样，根据一元线性函数的表达式$  y(x, w) = w_0 + w_1x $\n",
    "那么，转换成矩阵形式为：\n",
    "\n",
    "![image](equation.svg)\n",
    "\n",
    "即：\n",
    "\n",
    "![image](equation--1-.svg)\n",
    "\n",
    "然后，平方损失函数转换为矩阵形式：\n",
    "\n",
    "![image3](equation--4-.svg)\n",
    "\n",
    "此时，对矩阵求偏导数得到：\n",
    "\n",
    "![image3](equation--5-.svg)\n",
    "\n",
    "当矩阵$ X^TX $满秩（不满秩无法通过 OLS 方法求解）时，$ (X^TX)^{-1}X^TX=E $，且$ EW=W $。所以，$ (X^TX)^{-1}X^TXW=(X^TX)^{-1}X^Ty $。最终得到：\n",
    "\n",
    "![image3](equation--6-.svg)\n",
    "\n",
    "以上，我们完成了线性回归矩阵求解推导。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[44.25604341],\n",
       "        [ 0.7175629 ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#最小二乘法矩阵求解实现\n",
    "#使用矩阵方法求解时的代码实现过程会更加简单，因为最终参数的求解就只有一个公式\n",
    "\n",
    "def ols_matrix(x, y):\n",
    "    w = (x.T * x).I * x.T * y\n",
    "    return w\n",
    "#由于是矩阵运算，所以需要对原数据进行形状变换，并添加截距项系数 1，否则会因为形状问题无法计算\n",
    "x = np.array(x).reshape(len(x), 1)\n",
    "x = np.concatenate((np.ones_like(x), x), axis=1)\n",
    "x = np.matrix(x)\n",
    "\n",
    "y = np.array(y).reshape(len(y), 1)\n",
    "y = np.matrix(y)\n",
    "\n",
    "ols_matrix(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面我们使用最小二乘法来求解线性方程系数，得到的结果是精确的解析解。尝试将拟合的直线绘制成图，查看拟合效果。四舍五入保留参数的 3 为小数。\n",
    "$ y=44.256+0.718∗x $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1d129e250c8>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAf8ElEQVR4nO3deXiU1fn/8fctuKBWI5vFAMW2LAoqYKrgWkopy9dCwFpBWkCtWIv+RBEhIiC0gIgVpS4VgaI/Fbci4AKIuOCGrCJhE75SIIlsIigQMZDz/eNMbIQJWWaSZ57J53VduSbzZCZzP9dwfTi5nzPnmHMOERFJLscEXYCIiMSfwl1EJAkp3EVEkpDCXUQkCSncRUSSUNWgCwCoWbOma9CgQdBliIiEytKlS3c652pF+1lChHuDBg1YsmRJ0GWIiISKmW0q6mdqy4iIJKFiw93MppjZdjPLLHTsr2b2qZl9YmZvmNkZkeNmZhPMbEPk5y3Ls3gREYmuJCP3qUCHw46Nc86d65xrDrwKDIsc7wg0jHz1BR6LU50iIlIKxYa7c24BsOuwY18XunsSULCGQRfgKectBFLMrE68ihURkZIp8wVVMxsF9AL2AG0ih1OBLYUelhU59kWU5/fFj+6pX79+WcsQEZEoynxB1Tk3xDlXD3gGuDly2KI9tIjnT3TOpTnn0mrVijqTR0SkWDOWZ3PxvW9x5uDXuPjet5ixPDvokhJCPGbLPAtcGfk+C6hX6Gd1gZw4vIaIyBFmLM8mY/pKsnfn4oDs3blkTF+pgKeM4W5mDQvd7QysjXw/C+gVmTXTCtjjnDuiJSMiEg/j5q4jN+/QD47l5h1i3Nx1AVWUOIrtuZvZNOCXQE0zywKGA53MrDGQD2wC/hx5+OtAJ2ADsB+4thxqFhEBIGd3bqmOVybFhrtzrkeUw5OLeKwD+sValIhISZyRUo3sKEF+Rkq1AKpJLPqEqoiE1sD2jal2bJUfHKt2bBUGtm8cUEWJIyHWlhERKYv0FqmA773n7M7ljJRqDGzf+PvjlZnCXURCLb1FqsI8CrVlRESSkMJdRCQJKdxFRJKQwl1EJAkp3EVEkpDCXUQkCSncRUSSkMJdRCQJKdxFRJKQwl1EJAkp3EVEkpDCXUQkCSncRUSSkMJdRCQJKdxFRJKQwl1EJAkp3EVEkpDCXUQkKJs2wa5d5fKrFe4iIhVt714YOhSaNIERI8rlJbSHqohIRcnPh6efhowMyMmBa66BO+4ol5fSyF1EpCJ8+CG0agW9e0Pduv7+M89AvXrl8nIKdxGR8rR5sx+hX3wxZGfDU0/BRx9B69bl+rJqy4iIlId9+2DsWBg3zt8fOhQGDYKTTqqQl1e4i4jEU34+PPssDB7sR+rdu/uQr1+/QstQW0ZEJF4WLoSLLoI//hHq1IH334dp0yo82EHhLiISuy1boGdP30ffvBmmToWPP/Z99oCoLSMiUlb79/ue+tixvh0zZIhvx5x8ctCVKdxFRErNOd9uGTQIsrLg97/3Ad+gQdCVfU9tGRGR0li0yPfVe/aE2rVhwQJ4/vmECnYoQbib2RQz225mmYWOjTOztWb2qZm9bGYphX6WYWYbzGydmbUvr8JFRCpUdjb06gUXXgj/+Q9MmQKLF8OllwZdWVQlGblPBTocdmwe0Mw5dy7wGZABYGZnA92BppHnPGpmVeJWrYhIRcvNhb/+FRo1ghde8EsHfPYZXHstHJO4zY9iK3POLQB2HXbsDefcwcjdhUDdyPddgOeccweccxuBDcAFcaxXRKRiOAfPPQeNG8OwYdCxI6xZA6NHw49+FHR1xYrHfzvXAbMj36cCWwr9LCty7Ahm1tfMlpjZkh07dsShDBGROFm8GC65BHr0gBo14J134KWX4Mwzg66sxGIKdzMbAhwEnik4FOVhLtpznXMTnXNpzrm0WrVqxVKGiEh85ORAnz5wwQWwYQNMmgRLlsDllwddWamVeSqkmfUGrgDaOucKAjwLKLzEWV0gp+zliYhUgNxceOABGDMG8vLgzjv9nPVTTgm6sjIr08jdzDoAg4DOzrn9hX40C+huZseb2ZlAQ2BR7GWKiJQD5/xF0rPOgrvvht/8Blav9nPWQxzsUIKRu5lNA34J1DSzLGA4fnbM8cA8MwNY6Jz7s3NulZm9AKzGt2v6OecOlVfxIiJltmwZ9O8P770H554Lb70FbdoEXVXcFBvuzrkeUQ5PPsrjRwGjYilKRKTcbN0Kd93l13+pWRMmToTrroMqyTVrW8sPiEjl8O23MH68n8p44IDf3m7IEDj11KArKxcKdxFJbs7B9OkwcCBs3AhdusD998PPfx50ZeUqcT9eJSISq+XLfR/9d7/zKzW++SbMmJH0wQ4KdxFJRtu2wQ03wPnnw6pV8M9/+guobdsGXVmFUVtGRJLHgQPw0EPwt7/5ueu33eb3Lk1JKf65SUbhLiLh55xvt9xxB3z+Ofz2t76v3qhR0JUFRm0ZEQm3FSt8u6VbN6hWDd54A2bNqtTBDgp3EQmr7duhb19o0QI+/RQeeQQ++QTatQu6soSgtoyIhMuBA/CPf/g11vfvh1tv9UvynnZa0JUlFIW7iISDc77dMmAA/O//wv/8j++rN2kSdGUJSW0ZEUl8K1f6dkt6Ohx3HMyZA6++qmA/CoW7iCSuHTvgppugeXP/gaR//MNfQG2v7ZmLo7aMiCSe776Dhx+GkSNh7164+WYYPhyqVw+6stBQuItI4nDOt1sGDID166FDB7+JxllnBV1Z6KgtIyKJITPTt1s6d/bL777+OsyerWAvI4W7iARr507o1w/OO89vTP3QQ37eeseOQVcWamrLiEgw8vL8B49GjIBvvvEXTkeMgBo1gq4sKSjcRaRiOedbLgMGwLp1ft/SBx6Apk2DriypqC0jIhVn9Wrfbrniiv9ePJ0zR8FeDhTuIlL+vvwSbrnFb0S9cKEfqa9c6T9lahZ0dUlJbRkRKT95efDYY3DPPbBnD9x4o5+7XrNm0JUlPYW7iJSPOXPg9tthzRr49a/95tTNmgVdVaWhtoyIxNfatdCpk++t5+XBzJl+jXUFe4VSuItIfOza5ZffPecc+OADv2LjqlX+Q0nqq1c4tWVEJDYHD8Ljj/s11Xfv9htTjxwJtWsHXVmlppG7iJTdG2/4T5befLO/Xb4c/vlPBXsCULiLSOmtW+c3oW7f3u+M9PLLMH++n+ooCUHhLiIl99VXfgZMs2bw7rtw332+r56err56glHPXUSKd/AgPPEEDB3qL5z+6U9+D9PTTw+6MimCRu4icnRvvgktWsBf/uJH7MuWwcSJCvYEp3AXkejWr/fTGNu1g3374N//hrff9lveScIrNtzNbIqZbTezzELHrjKzVWaWb2Zphz0+w8w2mNk6M9NGhyJhs3u3X7GxaVMf5vfe6xf86tZNffUQKcnIfSrQ4bBjmUA3YEHhg2Z2NtAdaBp5zqNmViX2MkWk3B065OerN2zolwro1cuP3gcNghNOCLo6KaViw905twDYddixNc65dVEe3gV4zjl3wDm3EdgAXBCXSkWk/Lz1lu+r//nPflu7JUtg0iT48Y+DrkzKKN4991RgS6H7WZFjRzCzvma2xMyW7NixI85liEiJbNgAXbtC27Z+N6QXX/RTHFu2DLoyiVG8p0JGa8i5aA90zk0EJgKkpaVFfYyIeDOWZzNu7jpydudyRko1BrZvTHqLqOOmktmzB0aNggcfhOOOg9Gj4bbb1H5JIvEO9yygXqH7dYGcOL+GSKUyY3k2GdNXkpt3CIDs3blkTF8JUPqAP3QIpkyBu++GHTugTx8f8nXqxLlqCVq82zKzgO5mdryZnQk0BBbF+TVEKpVxc9d9H+wFcvMOMW5utMteR/HOO3D++dC3r79ounixD3oFe1IqyVTIacBHQGMzyzKz682sq5llAa2B18xsLoBzbhXwArAamAP0c84dKup3i0jxcnbnlur4ET7/HK68Etq08csHPP88vPeeD3pJWsW2ZZxzPYr40ctFPH4UMCqWokTkv85IqUZ2lCA/I6Xa0Z/49de+lz5+PFSt6pcLGDAAqhXzPEkK+oSqSIIb2L4x1Y794cdFqh1bhYHtG0d/wqFDMHkyNGoEY8dC9+5+vvrddyvYKxEtHCaS4AoumpZotsyCBdC/v19X/aKL4JVX4Be/qOCKJREo3EVCIL1F6tFnxmzcCHfeCS+9BPXqwbPP+hG7lguotBTuImH2zTcwZgw88ABUqQIjRsAdd8CJJwZdmQRM4S4SRvn58NRTkJEBW7fCH/7gQ75u3aArkwShcBcJm/ff9331pUvhwgv9FnetWgVdlSQYzZYRCYtNm+Dqq+HSS/1o/emn4cMPFewSlUbuIolu714/pfH++/0F0uHDYeBAOOmkoCuTBKZwF0lU+fl+dJ6RATk5cM01fuOMevWKf65UemrLiCSignZL796QmurvP/OMgl1KTOEukkg2b/Yj9IsvhuxsPyNm4UJo3TroyiRk1JYRSQT79sF998G4ceAcDB3qP5R08slBVyYhpXAXCVJ+vv806eDBfqR+9dX+4ulPfhJ0ZRJyasuIBGXhQr/+yx//6Pcqff99eO45BbvEhcJdpKJt2QI9e/o++ubNMHUqLFrk++wicaK2jEhF2b/f99THjvXtmLvu8tMcE7yvHvf9W6VCKNxFyptzMG0aDBoEWVlw1VX+4mmDBkFXVqy47t8qFUptGZHytGiR76v37Am1asG778ILL4Qi2CGO+7dKhVO4i5SH7Gzo1csv7LVxo98ZafFiuOyyoCsrlZj3b5XAKNxF4ik31+9V2qiR34h68GC/xd111/n11kOmqH1ai92/VQKncBeJB+f8NMbGjWHYMOjYEdau9Wus/+hHQVdXZqXev1UShsJdJFaLF8Mll0CPHlCjBrzzjt/u7swzg64sZuktUhnT7RxSU6phQGpKNcZ0O0cXU0NAs2VEyionx09nfPJJqF0bJk2CPn1C2X45mmL3b5WEpHAXKa3cXL9n6ZgxkJfn14AZMgROOSXoykS+p3AXKSnn4MUXfZhv2gRdu/oPJf3sZ0FXJnIE9dxFSmLZMrj8cr+w16mnwltvwfTpCnZJWAp3kaPZutVPY0xL87NfHn/cB32bNkFXJnJUasuIRPPttzB+PIweDQcOwIABcPfdftQuEgIKdyl3ibTwVLG1OOfbLQMH+k+Wduni++oNGwZSr0hZKdylXCXSwlPF1rJ8Odx2m1//pVkzePNNaNu2QmsUiRf13KVcJdLCU0XVMumlj+CGG+D882HVKnjsMR/0CnYJsWJH7mY2BbgC2O6caxY5Vh14HmgA/Af4vXPuKzMz4CGgE7Af6OOcW1Y+pUsYJNLCU4e/5nEH87h26Uxu/vB5yM/zo/ahQyElpcJrE4m3kozcpwIdDjs2GJjvnGsIzI/cB+gINIx89QUei0+ZElaJtPDU96/pHO0/+5B5k28i452prPhZcz9i//vfFeySNIoNd+fcAmDXYYe7AE9Gvn8SSC90/CnnLQRSzKxOvIqV8EmkhacGtm9M812bePa5ITz+8mi+rXoc1/f4GzufedGv4iiSRMp6QfV059wXAM65L8ysduR4KrCl0OOyIse+KHuJEmYFF00Dny2zfTvp/xxJl8mT2HPCyQxtdxPvXt6V2zudrXVTJCnFe7aMRTnmoj7QrC++dUP9+vXjXIYkkkAXnvruO5gwwa+xvn8/dsstpAwfzl9POy2YekQqSFlny2wraLdEbrdHjmcB9Qo9ri6QE+0XOOcmOufSnHNptWrVKmMZIkVwDmbOhKZN/Zz1Sy6BlSvhwQdBwS6VQFnDfRbQO/J9b2BmoeO9zGsF7Clo34hUmJUroV07SE+HY4+F2bPhtdegSZOgKxOpMMWGu5lNAz4CGptZlpldD9wLtDOz9UC7yH2A14HPgQ3AE8BfyqVqkWh27ICbboLmzf36LxMmwIoV0OHwyV4iya/YnrtzrkcRPzriEx7OOQf0i7UokVL57jt4+GEYORL27oV+/eCee6B69aArEwmMlh+Q8HIOXn3VL+q1fr0foT/wAJx1VtCVxV0irc8j4aDlByScMjOhfXvo3BmOOcb31GfPTtpgz5i+kuzduTj+uybOjOXZQZcmCUzhLuGyc6dvu5x3nt+Y+sEH/QXUTp2CrqzcJNL6PBIeastIOOTlwSOPwIgR8M03/sLpPfdAzZpBV1buEml9HgkPjdwlsTnnWy7nnOMX9vrFL/wMmIcfrhTBDom1Po+Eh8JdEtfq1dCxI1xxBeTnwyuvwNy5/oNJlUgirc8j4aFwl8Tz5Zdwyy1w7rmwcKGfAZOZ6UPeoq1wkdzSW6Qypts5pKZUw4DUlGqM6XaOZsvIUannLokjL89vlHHPPbBnD9x4o++xa3mKYNfnkVBSuEtimDMHbr8d1qzxOyCNH+/77CJSJmrLSLDWrvXTGDt29CP3mTNh3jwFu0iMFO4SjF274NZbfYh/8AHcf7/fDalz50rZVxeJN7VlpGIdPAiPPw7DhsHu3X5j6pEjoXbt4p8rIiWmcJeK88Ybfq766tXQpo3/dOm555b612idFZHiqS0j5W/dOvjtb/1aMAcOwMsvw/z5ZQ52rbMiUjyFu5Sfr77yM2CaNYN334X77vN99fT0MvfVtc6KSMmoLSPxd/AgPPEEDB3qL5z+6U9+D9PTT4/5V2udFZGS0chd4uvNN6FFC/jLX/yIfdkymDgxLsEOWmdFpKQU7hIf69dDly5+79J9++Df/4a33/Zb3sWR1lkRKRm1ZSQ2e/b4lsuECXD88TBmDPTvDyecUC4vVzArRrNlRI5O4S5lc+gQTJrk++o7d8K118KoUfDjH5f7S2udFZHiKdyl9N56y4/OV66ESy/189Vbtoz512r+ukj8qOcuJbdhA3Tt6hf2+vpreOEFP8UxTsGu+esi8aNwl+Lt2QN33glnn+0X9Ro1yi/4ddVVcVsHRvPXReJLbRkp2qFDMGUK3H03bN8OffrA6NFQp07cX0rz10XiSyN3ie6dd+D886FvX2jYEBYvhn/9q1yCHTR/XSTeFO7yQ59/Dlde6Rf2+uoreO45eO89SEsr15fV/HWR+FJbRryvv/Ytl/HjoWpVP3d9wACoVjEjZ81fF4kvhXtld+gQTJ0KQ4bAtm3Qq5cP+dSKD1XNXxeJH4V7MZJ67vWCBX6++vLl0Lo1zJoFF1wQdFUiEgfquR9F0s693rjRT2O8/HLYsQOefdZvdadgF0kaCvejSLq51998A3fdBWedBa+9BiNG+I00evTQvqUiSUZtmaNImrnX+fnw1FOQkQFbt0LPnnDvvVC3btCViUg5iWnkbma3mlmmma0ys/6RY9XNbJ6ZrY/cnhafUiteUsy9fv9932659lr4yU/go4/g6acV7CJJrszhbmbNgBuAC4DzgCvMrCEwGJjvnGsIzI/cD6VQz73etAmuvtov7LV1qw/0Dz+EVq2CrkxEKkAsI/ezgIXOuf3OuYPAu0BXoAvwZOQxTwLpsZUYnPQWqYzpdg6pKdUwIDWlGmO6nZPYs2X27vXL8DZp4me/DBvm++o9e8IxusQiUlmYc65sTzQ7C5gJtAZy8aP0JcAfnXMphR73lXPuiNaMmfUF+gLUr1///E2bNpWpDonIz/ej84wMyMnxF0nvvRfq1w+6MhEpJ2a21DkX9ePjZR7KOefWAGOBecAcYAVwsBTPn+icS3POpdWqVausZQj8t93Su7f/8NEHH/jpjQp2kUorpr/TnXOTnXMtnXOXAbuA9cA2M6sDELndHnuZEtXmzXDNNXDxxZCVBU8+CQsXwkUXBV2ZiAQs1tkytSO39YFuwDRgFtA78pDe+NaNxNO+fTB8uO+rv/yyX5L3s8/80gHqq4sIsc9z/7eZ1QDygH7Oua/M7F7gBTO7HtgMXBVrkRKRn+/bLYMHQ3a2nw0zdqyf4igiUkhM4e6cuzTKsS+BtrH8Xoli4UK/DszHH/t11p9/3rdjRESi0N/wiS4rC/7wB7+w1+bNfgXHRYsU7CJyVFp+IFHt3w/jxvm2S36+XxMmIwNOPjnoykQkBBTuicY5mDYNBg3yo/arroL77oMGDYKuTERCROGeSBYtgltv9f31Fi3gmWfgssuCrgpI8nXtRZKQeu6JIDvbT2O88EK/1vrkyX5D6gQK9qRc114kiSncg5Sb6/cqbdTIz34ZPBjWr4frroMqVYp/fgVJunXtRSoBtWWC4JwP8zvvhC1b4MorfV/9pz8NurKokmZde5FKRCP3irZ4MVxyiV/Yq3p1ePtteOmlhA12SJJ17UUqGYV7RcnJgT59/MYZGzbAE0/A0qXwy18GXVmxQr2uvUglpbZMecvNhQcegDFjIC/Pt2KGDIFTTgm6shIrmBWj2TIi4aFwLy/OwYsv+jDftAm6dvUfSvrZz4KurEzSW6QqzEVCRG2Z8rBsGVx+uV/Y69RTYf58mD49tMEuIuGjcI+nrVv9NMa0NFi7Fh5/3Af9r34VdGUiUsmoLRMP334L48fD6NFw4ADcfrvfx/TUU4OuTEQqKYV7LJzz7ZaBA/0nSzt3hvvvh4YNg65MRCo5tWXKavlyaNMGfvc7OOkkmDcPZs5UsItIQlC4l9a2bXDDDX7DjMxMePRRH/S//nXQlYmIfE9tmZI6cAAeegj+9jc/d71/fxg2DFJSgq5MROQICvfiOAczZsAdd8Dnn8MVV8Df/+4X+xIRSVBqyxzNihXQti106wYnnABz58IrryjYRSThKdyj2b4dbrwRWrb0Af/ww/72N78JujIRkRJRW6aw776DCRP8Guv798Mtt/i+evXqQVcmIlIqCnfwffVZs3xffcMG6NTJ99WbNAm6MhGRMlFbZuVKaNcO0tOhalWYPRtee03BLiKhVnnDfccOuOkmaN7cr/8yYQJ8+il06BB0ZSIiMat8bZnvvvMXSEeOhL17oV8/GD4catQIujIRkbipPOHuHLz6KgwY4Dehbt/eb6Jx9tlBVyYiEneVoy2TmenDvHNnOOYY31OfPVvBLiJJK7nDfedO33Y57zy/MfWDD/oLqJ06gVnQ1YmIlJvkbMvk5cEjj8CIEfDNN/7C6T33QM2aQVcmIlIhkivcnYPXX/d99XXr/BTH8eOhadOgKxMRqVAxtWXM7DYzW2VmmWY2zcxOMLMzzexjM1tvZs+b2XHxKvaoVq+Gjh39wl75+X4NmLlzFewiUimVOdzNLBX4f0Cac64ZUAXoDowFxjvnGgJfAdfHo9AiffmlXybg3HNh4UI/AyYz04e8+uoiUknFekG1KlDNzKoCJwJfAL8CXor8/EkgPcbXKNrrr/udjx59FPr29VMcb7sNjquYPxZERBJVmcPdOZcN3A9sxof6HmApsNs5dzDysCwgNdYii9SoEbRqBZ984gO+Vq1yeykRkTAp8wVVMzsN6AKcCewGXgQ6RnmoK+L5fYG+APXr1y9bET//uR+9i4jID8TSlvk1sNE5t8M5lwdMBy4CUiJtGoC6QE60JzvnJjrn0pxzabU04hYRiatYwn0z0MrMTjQzA9oCq4G3gd9FHtMbmBlbiSIiUlqx9Nw/xl84XQasjPyuicAg4HYz2wDUACbHoU4RESmFmD7E5JwbDgw/7PDnwAWx/F4REYlNcq8tIyJSSSncRUSSkMJdRCQJKdxFRJKQORf1M0YVW4TZDmBTGZ9eE9gZx3KCpHNJTMlyLslyHqBzKfAT51zUDwolRLjHwsyWOOfSgq4jHnQuiSlZziVZzgN0LiWhtoyISBJSuIuIJKFkCPeJQRcQRzqXxJQs55Is5wE6l2KFvucuIiJHSoaRu4iIHEbhLiKShEIX7mb2HzNbaWafmNmSyLHqZjYvsin3vMhGIgnPzFLM7CUzW2tma8ysddjOxcwaR96Lgq+vzax/2M6jQEJt+h4jM7s1ch6rzKx/5Fgo3hczm2Jm280ss9CxqLWbN8HMNpjZp2bWMrjKf6iI87gq8p7km1naYY/PiJzHOjNrH8trhy7cI9o455oXmhs6GJgf2ZR7fuR+GDwEzHHONQHOA9YQsnNxzq2LvBfNgfOB/cDLhOw8IIE2fY8DM2sG3IBfofU84Aoza0h43pepQIfDjhVVe0egYeSrL/BYBdVYElM58jwygW7AgsIHzexs/L+3ppHnPGpmVcr8ys65UH0B/wFqHnZsHVAn8n0dYF3QdZbgPE4BNhK5qB3mcylU+2+AD8J6Hvj9frcA1fHLYb8KtMd/erBq5DGtgblB11qCc7kKmFTo/lDgzjC9L0ADILPQ/ai1A48DPaI9LhG+Dj+PQsffwQ8kCu5nABmF7s8FWpf1dcM4cnfAG2a2NLIPK8DpzrkvACK3tQOrruR+CuwA/mVmy81skpmdRDjPpUB3YFrk+9Cdh0uETd/jJxO4zMxqmNmJQCegHiF8XwopqvaC/5QLhOU9OlxczyOM4X6xc64l/k+xfmZ2WdAFlVFVoCXwmHOuBbCPxP0TuViRPnRn/EbpoXTYpu9nACdRik3fE4lzbg2+nTQPmAOsAA4e9UnhZVGOJfx7FEVczyN04e6cy4ncbsf3di8AtplZHYDI7fbgKiyxLCDL+e0KwW9Z2JJwngv4EFzmnNsWuR/G84hp0/dE45yb7Jxr6Zy7DNgFrCec70uBomrPwv9VUiA079Fh4noeoQp3MzvJzH5U8D2+x5sJzMJvxg0h2ZTbObcV2GJmjSOHCjYYD925RPTgvy0ZCOd5JNWm72ZWO3JbH38BbxrhfF8KFFX7LKBXZNZMK2BPQfsmZGYB3c3seDM7E3+BeFGZf1vQFxtKeWHip/g/L1cAq4AhkeM18FfP10duqwddawnPpzmwBPgUmAGcFsZzAU4EvgROLXQsdOcRqXsEsBY/aPj/wPGRf3eLgA34ttPxQddZwnN5D/+f0wqgbZjeF/x/RF8AefgR7fVF1Y5vZzwC/C+wkkIXKYP+KuI8uka+PwBso9AFemBI5DzWAR1jeW0tPyAikoRC1ZYREZGSUbiLiCQhhbuISBJSuIuIJCGFu4hIElK4i4gkIYW7iEgS+j9BkwDfDR/0xAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "x = [55, 71, 68, 87, 101, 87, 75, 78, 93, 73]\n",
    "y = [91, 101, 87, 109, 129, 98, 95, 101, 104, 93]\n",
    "\n",
    "def ols_algebra(x, y):\n",
    "    n = len(x)\n",
    "    w1 = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x*x) - sum(x)*sum(x))\n",
    "    w0 = (sum(x*x)*sum(y) - sum(x)*sum(x*y))/(n*sum(x*x)-sum(x)*sum(x))\n",
    "    return w0, w1\n",
    "\n",
    "ols_algebra(np.array(x), np.array(y))\n",
    "\n",
    "plt.scatter(x, y)\n",
    "plt.plot(np.array([50, 110]), np.array([50, 110]) * 0.718 + 44.256, 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. 最小二乘法的局限性和适用场景　　\n",
    "![image5](5.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 梯度下降法求解数学推导\n",
    "介绍完最小二乘法，下面来说一说梯度下降法的求解过程，关于梯度下降法的原理和过程这里不再解释。梯度下降法作为迭代优化方法，需要提前定义好损失函数。这里选择均方误差 MSE 作为损失函数。\n",
    "\n",
    "![image](equation--13-.svg)\n",
    "\n",
    "其中，$ y_{i} $表示真实值，$ \\hat y_{i} $表示近似值，nn 则表示值的个数。近似值则通过当前已知计算，那么：\n",
    "\n",
    "![image](equation--14-.svg)\n",
    "\n",
    "此时，我们求解参数和截距项对应的梯度：\n",
    "\n",
    "![image](equation--15-.svg)\n",
    "\n",
    "那么，根据梯度就可以执行迭代更新，从而找到最优参数：\n",
    "\n",
    "![image](equation--16-.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.264673601681181, 0.017892580143733863)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#梯度下降法求解代码实现\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "x = [55, 71, 68, 87, 101, 87, 75, 78, 93, 73]\n",
    "y = [91, 101, 87, 109, 129, 98, 95, 101, 104, 93]\n",
    "\n",
    "def gradient_descent(x, y, lr, num_iter):\n",
    "    w1 = 0  # 初始参数为 0\n",
    "    w0 = 0  # 初始参数为 0\n",
    "\n",
    "    for i in range(num_iter):  # 梯度下降迭代\n",
    "        # 计算近似值\n",
    "        y_hat = (w1 * x) + w0\n",
    "        # 计算参数对应梯度\n",
    "        w1_gradient = -(2/len(x)) * sum(x * (y - y_hat))\n",
    "        w0_gradient = -(2/len(x)) * sum(y - y_hat)\n",
    "        # 根据梯度更新参数\n",
    "        w1 -= lr * w1_gradient\n",
    "        w0 -= lr * w0_gradient\n",
    "\n",
    "    return w1, w0\n",
    "\n",
    "w1, w0 = gradient_descent(np.array(x), np.array(y), lr=0.00001, num_iter=100)\n",
    "w1, w0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1d129ec7688>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9fX/8depC4JWUYjKooIVKTtCRK2VWpGCS4G69If9WmlFsQUX3KXU8rWtguJSsIoiKGgVi+wIiIALiggEEAgggqgQQImyqSAEOL8/PsPXEBMSMpPcmcn7+XjMY2bu3Mmc68STD5977vmYuyMiIunlR1EHICIiiafkLiKShpTcRUTSkJK7iEgaUnIXEUlDh0YdAED16tW9Tp06UYchIpJS5s+f/6W7ZxT2WlIk9zp16pCVlRV1GCIiKcXMPivqNU3LiIikISV3EZE0pOQuIpKGik3uZvasmW00s+xCXrvDzNzMqseem5kNNLNVZrbYzFqURdAiInJgJRm5DwPaF9xoZicBbYE1+TZfBNSL3boBg+IPUUREDlaxyd3dZwKbCnnpMeAuIH/nsY7A8x68D1Q1sxoJiVREREqsVHPuZtYBWOfuiwq8VAtYm+95TmxbYT+jm5llmVlWbm5uacIQEZEiHHRyN7MqQG/gb4W9XMi2QnsKu/tgd89098yMjEJr8EVE0ldeHvTrB/PmlcmPL83I/SdAXWCRmX0K1AYWmNmJhJH6Sfn2rQ2sjzdIEZG0snAhnHUW9OoFo0eXyUccdHJ39yXufry713H3OoSE3sLdPwcmANfEqmbOBra6+4bEhiwikqK++w5694Yzz4T162HUqDB6LwMlKYUcAcwG6ptZjpl1PcDuk4HVwCrgGaB7QqIUEUl1s2ZB8+bwwANwzTWwfDlcfnmZfVyxvWXc/apiXq+T77EDPeIPS0QkTXz9NfzlL/DEE3DyyTB1KvzqV2X+sbpCVUSkrEydCo0bh8R+002QnV0uiR2U3EVEEm/TJvjDH6B9e6hSBd55BwYMgKOOKrcQlNxFRBJp9Gho2BD+859w8nThQjj33HIPIyn6uYuIpLwNG+DGG2HMGGjRAl57LZxAjYhG7iIi8XCH554Lo/VJk0Jp45w5kSZ20MhdRKT0Pv0UunWDadPgvPNgyBA4/fSoowI0chcROXh79sDAgaESZvbsUA3z1ltJk9hBI3cRkYOzfDlcdx28916ohnn66VC/nmQ0chcRKYm8PLj//jCX/uGH8MILMHlyUiZ20MhdRKR48+dD166waBH89rfw+ONw/PFRR3VAGrmLiBRlxw64557QwXHjRhg7Fv7736RP7KCRu4hI4WbODHPrK1eGUfvDD0PVqlFHVWIauYuI5LdtG/ToAb/4BezeDdOnhxLHFErsoOQuIvK9KVNCeeOgQdCzJyxZAm3aRB1VqWhaRkTkq6/g1ltDBUzDhqHM8eyzo44qLhq5i0jF5Q4jR0KDBjBiBPztb7BgQcondtDIXUQqqvXroXt3GD8eMjPD3HrTplFHlTAlWWbvWTPbaGbZ+bb9w8wWm9kHZva6mdWMbTczG2hmq2KvtyjL4EVEDpo7DB0apl+mToX+/UMLgTRK7FCyaZlhQPsC2/q7e1N3bw68Cvwttv0ioF7s1g0YlKA4RUTit3o1XHhhKHFs3jycML3jDjg0/SYxik3u7j4T2FRg27Z8T48EPPa4I/C8B+8DVc2sRqKCFREplT174F//giZNYN48eOopeOMNOO20qCMrM6X+c2Vm9wPXAFuBX8Y21wLW5tstJ7ZtQyHv70YY3XNykvZmEJHkN27hOvpPXcH6LTuoWbUyd7arT6czan2/w9Kl4SKkOXPgkktCYq9dO7qAy0mpq2Xcvbe7nwS8CNwY22yF7VrE+we7e6a7Z2ZkZJQ2DBGpwMYtXEevMUtYt2UHDqzbsoNeY5YwbuE62LUL/v53OOMM+PhjeOklmDixQiR2SEwp5EvA5bHHOcBJ+V6rDaxPwGeIiPxA/6kr2JG3Z79tO/L2MHFIrAKmTx+44gpYtgyuugqssPFneirVtIyZ1XP3lbGnHYAPY48nADea2cvAWcBWd//BlIyISCKs37Jjv+dH5H3Hre++xHXzxkHNGjBhAvz61xFFF61ik7uZjQDOB6qbWQ7QB7jYzOoDe4HPgD/Fdp8MXAysArYDfyyDmEVEAKhZtTLrYgn+7DWL6fva49TdvIFxrS6h0+svwjHHRBxhdIpN7u5+VSGbhxaxrwM94g1KRKQk7mxXn/tfep+e04fwPx+8xqdVa9Dl6n785rarK3RiB12hKiIprNO6hbR7/kYOz93I4FaX8fIlXbn51832r5apoJTcRST15ObCLbfAiBFUbtwYJo6nW6tWobZaADUOE5FU4h4afDVsCKNGwX33hSXwWrWKOrKko5G7iKSGnBz485/h1VdDMh86NPRel0Jp5C4iyW3vXhg8GBo1ghkz4NFHQ791JfYD0shdRJLXqlVw/fXw1ltwwQUhyf/kJ1FHlRI0cheR5LN7d1iQukmTsHjGM8+EfutK7CWmkbuIJJclS0Kjr3nzoEMHePJJqKXSxoOlkbuIJIedO0MvmBYt4NNP4eWXYdw4JfZS0shdRKI3Z04YrS9dCldfDY89BtWrRx1VStPIXUSi8+23cNttcM45sHVrKHN84QUl9gTQyF1EovHGG6ESZvXqUL/erx8cfXTUUaUNjdxFpHxt2RKSeps2cMgh8Pbb4aSpEntCKbmLSPkZPz60Dnj2WbjrLli0CFq3jjqqtKTkLiJlb+NG6NwZOnWCjIxwAvXBB6Fy5agjS1tK7iJSdtzhP/+BBg1g7Fj4xz8gKyssgSdlSidURaRsrF0Lf/oTTJ4MZ58dGn01bBh1VBVGsSN3M3vWzDaaWXa+bf3N7EMzW2xmY82sar7XepnZKjNbYWbtyipwEUlSe/fCoEEhkb/1FgwYAO++q8RezkoyLTMMaF9g2zSgsbs3BT4CegGYWUOgM9Ao9p4nzeyQhEUrIsnto4/g/POhe/cwWs/OhptvDlUxUq6KTe7uPhPYVGDb6+6+O/b0faB27HFH4GV33+nunxAWylYXfZF0t3s3PPQQNGsWesM8+yy8/jrUrRt1ZBVWIk6oXgtMiT2uBazN91pObNsPmFk3M8sys6zc3NwEhCEikVi0CM46C+6+Gy66CJYtgz/+EcyijqxCiyu5m1lvYDfw4r5Nhezmhb3X3Qe7e6a7Z2ZkZMQThohEYedOuPfeUPmSkwOvvAKjR0ONGlFHJsRRLWNmXYBLgTbuvi+B5wAn5dutNrC+9OGJSFJ67z247jpYvhy6dIFHHoFq1aKOSvIp1cjdzNoDdwMd3H17vpcmAJ3NrJKZ1QXqAXPjD1NEksI338Att8DPfx6afr32GgwbpsSehIoduZvZCOB8oLqZ5QB9CNUxlYBpFubV3nf3P7n7UjMbCSwjTNf0cPc9ZRW8iJSjadOgW7fQa/3GG+GBB+DHP446KilCscnd3a8qZPPQA+x/P3B/PEGJSBLZvBluvx2eew7q14d33gkjd0lqaj8gIkUbOzZcfPT889CrF3zwgRJ7ilD7ARH5oc8/h5tuglGjoHlzmDQpLH8nKUMjdxH5njsMHx5G6xMnhnn1uXOV2FOQRu4iEnz2GdxwA0ydCueeC0OGwE9/GnVUUkoauYtUdHv3wr//DY0ahQZfjz8OM2cqsac4jdxFKrIVK6BrV5g1C9q1g6efhlNOiToqSQCN3EUqorw86Ns3NPpatixciDRlihJ7GtHIXaSiWbgQrr02lDVecUWYhjnxxKijkgTTyF2kovjuu1CrfuaZodRx9OjQ7EuJPS1p5C5SEbz7bphb/+ij0I73kUfg2GOjjkrKkEbuIuns669DH5jzzoNdu8ICGs8+q8ReASi5i6SrqVOhcWN48smw1N2SJdC2bdRRSTlRchdJN5s2hR7r7dtDlSphSmbAADjqqKgjk3Kk5C6SLtxDL5gGDeCll6B371AZ87OfRR2ZREAnVEXSwYYN0KNH6OLYokWYkmnePOqoJEIauYukMvfQZ71hw3AR0oMPwpw5SuyikbtIyvrkk7Ay0vTpoRpmyBA4/fSoo5IkUezI3cyeNbONZpadb9uVZrbUzPaaWWaB/XuZ2SozW2Fm7coiaJEKbc8eGDgwVMK8/36ohnnrLSV22U9JpmWGAe0LbMsGLgNm5t9oZg2BzkCj2HueNLND4g9TRIDQB+a888Ii1b/4BSxdCn/+M/xIM6yyv2J/I9x9JrCpwLbl7r6ikN07Ai+7+053/wRYBbRKSKQiFVleHvzzn3DGGeEq0xdeCKsjnXxy1JFJkkr0n/tawNp8z3Ni237AzLqZWZaZZeXm5iY4DJE0Mn8+Wxs1g3vvZeKpZ3HJn55mXKNfglnUkUkSS3RyL+y3zQvb0d0Hu3umu2dmZGQkOAyRNLBjB9x9N96qFTs3fMH1l/2VmzrezdLdR9BrzBLGLVwXdYSSxBJdLZMDnJTveW1gfYI/QyT9zZwJ110HK1cy8cyL+evPrmHbEd9fYbojbw/9p66g0xmF/sNYJOEj9wlAZzOrZGZ1gXrA3AR/hkj62rYNuncPJ0t374bp07nlgu77JfZ91m/ZEUGAkipKUgo5ApgN1DezHDPrama/MbMc4BxgkplNBXD3pcBIYBnwGtDD3feUXfgiaWTy5LCO6VNPwa23hkZfbdpQs2rlQncvarsIlGBaxt2vKuKlsUXsfz9wfzxBiVQoX34JPXvCiy+GK01feQXOPvv/Xr6zXX16jVnCjrzvx0mVDzuEO9vVjyJaSRG6QlUkKu4wciTcdBNs3gx9+oSVkipV2m+3ffPq/aeuYP2WHdSsWpk729XXfLsckJK7SBTWrw8XH02YAJmZMGMGNGlS5O6dzqilZC4HRZe1iZQn99ADpmHDsCrSww/D7NkHTOwipaGRu0h5+fjj0OjrjTdCNcyQIXDaaVFHJWlKI3eRsrZnDzz6aBidz5sHTz8dErwSu5QhjdxFylJ2NnTtCnPnwqWXwqBBULt21FFJBaCRu0hZ2LUL7rsvrIq0enVY9m7CBCV2KTcauYsk2rx5cO21YdT+u9/Bv/4FKdw/adzCdSrDTEEauYskyvbtcMcd4QKkzZvDSP3FF1M+sfcas4R1W3bgwLotO9S0LEUouYskwptvQtOm8MgjcP31YRGNX/866qji1n/qiv2ujIXvm5ZJclNyF4nH1q1www1wwQXh+RtvhN4wxxwTbVwJUlRzMjUtS35K7iKlNXFiuBhpyJAwHbN4Mfzyl1FHlVBqWpa6lNxFDlZubjhR2qEDVKsWFqnu3x+qVIk6soS7s119Kh+2/zLIalqWGlQtI1JS7jBiBNx8c+i7ft99cM89cPjhUUdWZtS0LHUpuYuURE5OaPT16qtw1lkwdGjovV4BqGlZalJylzKXTHXSBx3L3r3wzDNw551hZaRHHw0j90MOKfo9IklAyV3K1L466X3ldPvqpIFyT/AHHcvKlaGs8e23QzXMM8/AqaeWZ8gipVaSZfaeNbONZpadb9txZjbNzFbG7o+NbTczG2hmq8xssZm1KMvgJfklU510iWPZvTu04m3aFD74IFTDTJ+uxC4ppSTVMsOA9gW23QPMcPd6wIzYc4CLCIti1wO6AYMSE6akqmSqky5RLIsXwznnhGmYdu1g2bLQ+MusnKIUSYxik7u7zwQ2FdjcERgeezwc6JRv+/MevA9UNbMaiQpWUk8y1UkfMJadO8Mydy1bwmefwX//C2PHQs2a5RylSGKUts79BHffABC7Pz62vRawNt9+ObFtP2Bm3cwsy8yycnNzSxmGJLtkqpMuKpYHanwTujf+/e/QuTMsXw6//a1G65LSEn0RU2H/N3hhO7r7YHfPdPfMjBRurCQH1umMWvS9rAm1qlbGgFpVK9P3siaRVMsUjOUnVWDi6tH84o+dQt36pEnwwgvhwiSRFFfaapkvzKyGu2+ITbtsjG3PAU7Kt19tYH08AUrqS6Y66f+LZcaMUAnzySfQvTv07QtHHx11eEVKpnJSSQ2lHblPALrEHncBxufbfk2sauZsYOu+6RuRpLBlC1x3HVx4IRx6aChzfOKJpE/sarsrB6skpZAjgNlAfTPLMbOuQD+grZmtBNrGngNMBlYDq4BngO5lErVIaYwfHxp9DRsGd98NixZB69ZRR1WsZConldRR7LSMu19VxEttCtnXgR7xBiWSUF98Ea4qHTkSmjUL3Rxbtow6qhJLpnJSSR3qCinpyz2cIG3YEMaNg3/+MyyBl0KJHZKrnFRSh5K7pKc1a+CSS+Caa6B+/XClae/ecNhhUUd20JKpnFRSh3rLSHrZuzeshHT33eHxgAHQo0dKN/pS210pDSV3SR8ffRQqYd55B9q2hcGDoU6dqKNKiGQqJ5XUoOQuKadgzfddbX5CxzdeDu0DKleG556DLl10halUaEruklIKtu09ZsVSThtwA3y+Cn7zm1CzXkPtjER0QlVSyr6a70q7d3HHzOeZMLwnx3/9Jb3/pw+MGaPELhKjkbuklPVbdtAiZzkPTRnAaZtyGNW4Df+44Dq2Vf4x90cdnEgSUXKX1PHNNzw0cyiXzx7H+qOrc82V9zHz1FCzXks13yL7UXKXpHHA5livvw7dunHFmjW8mHkpfX/+e76tVAVQzbdIYZTcJSkUtb7pYdu2cMmwh0M/mPr1sZkzOerIulSduoLtqvkWKZKSuySFwppjtc5+h7MGPAXbt0KvXvC3v8ERR9CJ8l9cWyTVKLlLUsjfBCvjm83cN20QF3/0HkuPP5Xqb0+HM86IMDqR1KPkLkmhZtXKrNu8nSuyZ/DXN4ZQOW8nD/6iC5Pa/o6ZSuwiB03JXZJCnyZVqHLznfx89QLm1m7IPe1vZsOJp9D34kZRhyaSkpTci6HlzcrY3r3wxBP8qlcv8hwe6XATT/y0LTWOPZK++m8tUmpK7gdQVAUH6IReQnz4YWj0NWsWtGvHYU8/ze2nnMLtUcclkgbiaj9gZreYWbaZLTWznrFtx5nZNDNbGbs/NjGhlj8tb1ZG8vLggQfCqkjLlsHw4TBlCpxyStSRiaSNUid3M2sMXA+0ApoBl5pZPeAeYIa71wNmxJ6nJC1vVgYWLIBWrcLCGR06wPLlYUENdXAUSah4Ru4NgPfdfbu77wbeBn4DdASGx/YZDnSKL8ToaHmzBNqxI9Sqt2oFn38emny98gqccELUkYmkpXiSezbQ2syqmVkV4GLgJOAEd98AELs/vrA3m1k3M8sys6zc3Nw4wig7Wt4sQd59F5o3h379Qp/1ZctCe14RKTOlTu7uvhx4EJgGvAYsAnYfxPsHu3umu2dmZGSUNowy1emMWvS9rAm1qlbGCM2p+l7WRCdTS+rrr+HGG+G882DXLpg2DYYOhWNT9jSMSMqIq1rG3YcCQwHM7AEgB/jCzGq4+wYzqwFsjD/M6Gh5s1KaMgVuuAFycuCWW+Cf/4Sjjoo6KpEKI95qmeNj9ycDlwEjgAlAl9guXYDx8XyGpJivvgonSC++OCTzWbPgX/9SYhcpZ/HWuY82s2pAHtDD3TebWT9gpJl1BdYAV8YbpKQAdxg1KkzDbNoEf/1ruFWqFHVkIhVSvNMy5xWy7SugTTw/V5LPAa/U3bABuneHceOgZcvQe71Zs2gDFqngdIWqFKvIK3Xd6fTB63DbbbBzJzz0ENx6KxyqXyuRqOn/QilWYVfqVstdR63/1xtWLYDWreGZZ+D00yOKUEQKUnKXYuW/IvdHe/fQZcGr3DnzefbYj2DQIOjWDX4U17l5EUkwJXcpVs2qlVm3ZQenfbmGh6YMoMX6Fbx5aksGXnkHY//UOerwRKQQSu5SrLsuqMtnd/XhhndH8O3hVbjl0tt5vVkb+l7eNOrQRKQISu5yYFlZdLy+KyxezLSmv6RX665UqlVDvdZFkpySuxRuxw7o0wceeQROPBHGj6dthw60jTouESkRJXf5obffDotorFoF118fShyrVo06KhE5CCpxkO9t2wZ//jOcf35Y/m7GDBg8WIldJAUpuUswaRI0ahSS+W23weLFcMEFUUclIqWk5F7RffklXH01XHopHHMMvPdemGc/8sioIxOROCi5V1Tu8PLL0KABjBwZTp4uWABnnRV1ZCKSADqhWhGtWxcafU2YAGeeGRbQaNIk6qhEJIE0cq9I3EMPmIYNw6pIDz8Ms2crsYukIY3cK4qPPw5ljW++GaphnnkGTjst6qhEpIxo5J7u9uyBRx8No/P58+Hpp0OJoxK7SFrTyD2dZWdD164wd26ohhk0CGrXjjoqESkH8a6hequZLTWzbDMbYWZHmFldM5tjZivN7L9mdniigpUS2rUL7rsPWrSA1athxIhw8lSJXaTCKHVyN7NawM1Aprs3Bg4BOgMPAo+5ez1gM9A1EYFKCc2dG5a6+9//hSuvhOXLoXNnMIs6MhEpR/HOuR8KVDazQ4EqwAbgAmBU7PXhQKc4P0NKYvt2uP12OOcc2LwZJk6EF1+E6tWjjkxEIlDq5O7u64CHgTWEpL4VmA9scffdsd1ygEL7wppZNzPLMrOs3Nzc0oYhECpgmjQJJ06vvx6WLg1z7CJSYcUzLXMs0BGoC9QEjgQuKmRXL+z97j7Y3TPdPTMjI6O0YVRsW7eGJe4uuCBMu7z5Jjz1VGgjICIVWjzTMhcCn7h7rrvnAWOAnwFVY9M0ALWB9XHGKIWZODFcjDR0KNx5Z2j0df75UUclIkkinuS+BjjbzKqYmQFtgGXAm8AVsX26AOPjC1H2k5sLV10FHTpAtWowZ07ot16lStSRiUgSiWfOfQ7hxOkCYEnsZw0G7gZuM7NVQDVgaALiFPdwgrRBAxg9Gv7+d8jKgszMqCMTkSQU10VM7t4H6FNg82qgVTw/VwpYuzYsojFpUujaOHRo6L0uIlIEtR9IZnv3hhOkjRqFk6WPPQazZimxi0ix1H4gWa1cGcoa334b2rQJKySdemrUUYlIitDIPdns3g39+0PTpvDBB2EKZto0JXYROSgauSeTRYtCo6/586FjR3jySahZM+qoRCQFaeSeDHbuhHvvDZUva9eGZe/GjlViF5FS08g9arNnh9H68uXw+9+Hk6bVqkUdlYikOI3co/Ltt9CzJ5x7LnzzDUyeDM8/r8QuIgmhkXsUpk8PlTCffhoWqu7bF44+OuqoRCSNaORenrZsCVMwbdvCYYfBzJnwxBNK7CKScEru5WXcuNDoa/hwuOeeUBlz3nlRRyUiaUrTMmXtiy/gppvglVegWbPQzbFly6ijEpE0p5F7WXEPJ0gbNIDx4+H++2HePCV2ESkXGrmXhTVr4IYb4LXXwrJ3Q4eGJC8iUk40ck+kvXvDCdJGjeCdd2DgwHCvxC4i5Uwj90RZsQKuuw7efTdUwwweDHXqRB2ViFRQGrnHKy8P+vULJ0uzs+G552DqVCV2EYmURu7xWLgw1K0vXAiXXRamZE48MeqoRERKP3I3s/pm9kG+2zYz62lmx5nZNDNbGbs/NpEBJ4XvvoPeveHMM2H9ehg1Kix9p8QuIkkinjVUV7h7c3dvDrQEtgNjgXuAGe5eD5gRe54+Zs2C5s3hgQfg6qth2TK4/PKooxIR2U+i5tzbAB+7+2dAR2B4bPtwoFOCPiNa33wDN98crir97rtQ5jhsGBx3XNSRiYj8QKKSe2dgROzxCe6+ASB2f3xhbzCzbmaWZWZZubm5CQqjjEydCo0bw7//DTfeGE6ctmsXdVQiIkWKO7mb2eFAB+CVg3mfuw9290x3z8zIyIg3jLKxaRP84Q/Qvj0cccT3tetHHRV1ZCIiB5SIkftFwAJ3/yL2/AszqwEQu9+YgM8of6NHh0Zf//kP/OUvYT3Tc8+NOioRkRJJRHK/iu+nZAAmAF1ij7sA4xPwGeVnw4ZwgvSKK8Iyd1lZoS/MEUdEHZmISInFldzNrArQFhiTb3M/oK2ZrYy91i+ezyg37uEEacOGMGlSuDBpzpxQGSMikmLiuojJ3bcD1Qps+4pQPZM6Pv0UunWDadPg5z+HIUOgfv2ooxIRKbWK3X5gzx54/PFQCTN7drjC9O23ldhFJOVV3PYDy5eHRl/vvReqYZ56Ck45JeqoREQSouKN3PPywgnS5s3hww/DghqTJyuxi0haqVgj9wUL4Nprw/qlV14ZpmROOCHqqEREEq5ijNx37AiLUrdqFdY0HTMGRo5UYheRtJX+I/d33glz6x99FNrz9u8Px6Zfo0oRkfzSd+S+bRv06AGtW8OuXaHMccgQJXYRqRDSM7lPmRLKGwcNgp49Q6OvCy+MOioRkXKTXsn9q6/gmmvg4otDc69Zs+Cxx+DII6OOTESkXKVHcncPJ0gbNIARI+Dee8PSd+ecE3VkIiKRSP0TquvXh7n1ceOgZUuYPh2aNo06KhGRSKV2cp88GX73O9i5Ex56CG69FQ5N7UMSEUmE1M6Ep58epl4GDoR69aKORkQkaaR2cj/ttFAZIyIi+0mPE6oiIrIfJXcRkTSk5C4ikobiXWavqpmNMrMPzWy5mZ1jZseZ2TQzWxm71/X+IiLlLN6R+wDgNXf/KdAMWA7cA8xw93rAjNhzEREpR6VO7mZ2NNAaGArg7rvcfQvQERge22040CneIEVE5ODEM3I/FcgFnjOzhWY2xMyOBE5w9w0AsfvjC3uzmXUzsywzy8rNzY0jDBERKSie5H4o0AIY5O5nAN9yEFMw7j7Y3TPdPTMjIyOOMEREpCBz99K90exE4H13rxN7fh4huZ8GnO/uG8ysBvCWu9cv5mflAp+VKhCoDnxZyvcmGx1LckqXY0mX4wAdyz6nuHuho+NSX6Hq7p+b2Vozq+/uK4A2wLLYrQvQL3Y/vgQ/q9RDdzPLcvfM0r4/mehYklO6HEu6HAfoWEoi3vYDNwEvmtnhwGrgj4SpnpFm1hVYA1wZ52eIiMhBiiu5u/sHQGF/cdrE83NFRCQ+6XCF6uCoA0ggHUtySpdjSZfjAB1LsUp9QlVERJJXOozcRUSkACV3EZE0lHLJ3cw+NbMlZvaBmWXFtqVks7J0aLxmZvVj38W+2yEfYWEAAAOeSURBVDYz65lqx7GPmd1qZkvNLNvMRpjZEWZW18zmxI7lv7HqsKRnZrfEjmOpmfWMbUuJ78XMnjWzjWaWnW9bobFbMNDMVpnZYjNrEV3k+yviOK6MfSd7zSyzwP69YsexwszaxfPZKZfcY37p7s3z1YamarOylG+85u4rYt9Fc6AlsB0YS4odB4CZ1QJuBjLdvTFwCNAZeBB4LHYsm4Gu0UVZMmbWGLgeaEX43brUzOqROt/LMKB9gW1FxX4RUC926wYMKqcYS2IYPzyObOAyYGb+jWbWkPD71ij2nifN7JBSf7K7p9QN+BSoXmDbCqBG7HENYEXUcZbgOI4GPiF2UjuVjyVf7L8CZqXqcQC1gLXAcYQy4VeBdoSrBw+N7XMOMDXqWEtwLFcCQ/I9vxe4K5W+F6AOkJ3veaGxA08DVxW2XzLcCh5Hvu1vEQYS+573Anrlez4VOKe0n5uKI3cHXjez+WbWLbatRM3KkkxcjdeSVGdgROxxyh2Hu68DHiZcfLcB2ArMB7a4++7YbjmEPwLJLhtobWbVzKwKcDFwEin4veRTVOz7/ijvkyrfUUEJPY5UTO7nunsLwj/FephZ66gDKqW4Gq8lm9g8dAfglahjKa3YHG5HoC5QEziS8HtWUNLXD7v7csJ00jTgNWARsPuAb0pdVsi2pP+OCpHQ40i55O7u62P3Gwlzu62AL2JNyojdb4wuwhLLAXLcfU7s+ShCsk/FY4GQBBe4+xex56l4HBcCn7h7rrvnAWOAnwFVzWzf1dy1gfVRBXgw3H2ou7dw99bAJmAlqfm97FNU7DmEf5XskzLfUQEJPY6USu5mdqSZ/XjfY8IcbzYwgdCkDErYrCxq7v45sNbM9nXM3Nd4LeWOJeYqvp+SgdQ8jjXA2WZWxcyM77+TN4ErYvukyrFgZsfH7k8mnMAbQWp+L/sUFfsE4JpY1czZwNZ90zcpZgLQ2cwqmVldwgniuaX+aVGfbDjIExOnEv55uQhYCvSOba9GOHu+MnZ/XNSxlvB4mgNZwGJgHHBsKh4LUAX4Cjgm37aUO45Y3PcBHxIGDS8AlWK/d3OBVYRpp0pRx1nCY3mH8MdpEdAmlb4Xwh+iDUAeYUTbtajYCdMZTwAfA0vId5Iy6lsRx/Gb2OOdwBfkO0EP9I4dxwrgong+W+0HRETSUEpNy4iISMkouYuIpCEldxGRNKTkLiKShpTcRUTSkJK7iEgaUnIXEUlD/x+itQmjvpqHFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 因为是迭代方法，所以需要设置学习率和最大迭代次数。\n",
    "# 这里得到的结果与最小二乘法不太一致，这是不同方法带来的普遍现象\n",
    "# 同样绘制出拟合直线：\n",
    "\n",
    "plt.scatter(x, y)\n",
    "plt.plot(np.array([50, 110]), np.array([50, 110]) * w1 + w0, 'r')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "肉眼可见，最终的拟合情况不如最小二乘法。当然，如果随机初始参数以及调整合适的迭代次数及学习率，优化迭代结果可能会更理想一些。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k= 0.9004584204388926 b= 0.831055638876812\n",
      "cost：1\n",
      "求解的拟合直线为:\n",
      "y=0.9x+0.83\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 25311 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 21512 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 30452 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 32447 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 26679 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 26412 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 25968 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:211: RuntimeWarning: Glyph 25454 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 25311 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 21512 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 30452 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 32447 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 26679 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 26412 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 25968 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "C:\\Users\\PINE\\anaconda3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:180: RuntimeWarning: Glyph 25454 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3zV1f3H8dcRIoQIFRGjghZTFbAK0iaCiqtuJGIqCglLEHCgggtBakvVOuqouxaxiozQWr0/ja1Va92iJIQgVoJiohCQgGFJBiTk/P44QUWCkNxx7ng/Hw8e5A7u/XAfkHc+53uGsdYiIiIikbWX7wJEREQSkQJYRETEAwWwiIiIBwpgERERDxTAIiIiHiiARUREPGgZyTfbf//9bZcuXSL5liIiIt4sWLDga2ttx8Ye220AG2P+CvQH1lhrj2647x4gE9gKfA6MtNZu2N1rdenShYKCgqbULiIiErOMMV/u6rE9GYJ+GjjnB/e9Bhxtre0BfApMbnZ1IiIiCWi3AWytfRtY94P7XrXW1jXc/ADoHIbaRERE4lYoJmGNAl7e1YPGmLHGmAJjTMHatWtD8HYiIiKxL6gANsZMAeqA2bt6jrV2mrU23Vqb3rFjo9ehRUREEk6zZ0EbY0bgJmedbnWig4iISJM0K4CNMecANwGnWGurQluSiIhI/NvtELQxJheYB3Q1xpQZYy4FHgHaAq8ZY4qMMY+HuU4REZG4stsO2Fqb3cjdT4ahFhERkYShrShFREQ8UACLiIh4oAAWERHxQAEsIiIC8MwzcP31EXu7iJ6GJCIiEnW++QauvBJmzXK3Bw6E448P+9sqgEVEJHEtWACDB8OyZdCmDTz8MPTpE5G31hC0iIgknvp6uP9+1+kuWwY9ergwHjUKjIlICQpgERFJLGvWQP/+7npvbS1cfTV8+CF06xbRMjQELSIiieP112HoUFi9GvbbD/76VxgwwEsp6oBFRCT+1dbCzTfDmWe68D35ZFi0yFv4gjpgERGJd198AdnZ8MEHsNdeMHUqTJkCLVp4LUsBLCIi8evZZ2HMGNi4ETp3htmzXfcbBTQELSIi8aeqCi67DC6+2IXvBRe4IecoCV9QBywiIvFm8WK3tveTT6BVK7jvPrfRRoSWF+0pBbCIiMQHa+Hxx+G666Cmxi0rmjsXevb0XVmjNAQtIiKxb906t4XklVe68L30UigoiNrwBXXAIiIS6959F3JyYMUKaNcO/vIXNwQd5dQBi4hIbNq2DW67DU45xYXvccfBwoUxEb6gDlhERGLRypVuR6s333S3J06E22+HpCSvZTWFAlhERGLLSy/BJZdARQWkprpzfM86y3dVTaYhaBERiQ1btsCECZCZ6cL3rLPc2t4YDF9QBywiIrHg00/dtd2FC6FlS7jjDnea0V6x20cqgEVEJHpZ64aYx42DykpIS3NrezMyfFcWtNj90UFEROLbpk0wbJi73ltZ6Q5UWLgwLsIX1AGLiEg0KihwQ86ffw5t2sCjj8KIEVG3nWQw1AGLiEj0qK93ezefcIIL3549YcEC1wXHUfiCAlhERKLFmjVw3nlwww1QWwtXXeXO8O3WzXdlYaEhaBER8e8//3HXe1evhv32g6eegvPP911VWKkDFhERf2prYfJkt5Z39Wq3reSiRXEfvqAOWEREfCktdYcofPCBW887dSpMmQItWviuLCIUwCIiEnnPPgujR7ulRoccAnPmQN++vquKKA1Bi4hI5FRVwdixcPHFLnwvuACKihIufEEdsIiIV1W1VQSWBCjdUEpa+zSyumWRnJTsu6zwWLwYBg2CJUugVSu4/3644oq4W160pxTAIiKe5K/MJzM3k/LK8m/vS01JJS87j4xO8bHbE+C2k3z8cbj2WnegQvfubjvJHj18V+aVhqBFRDyorq3eKXwByivLyczNpLq22lNlIbZuHVx4IVx5pQvf0aMhPz/hwxcUwCIiXgSKAzuF73blleUEigMRrigM3n0Xjj0WAgFo1851vU88ASkpviuLCgpgEREPStaXBPV4VNu2DW67za3pXbECevd2E60GDfJdWVTRNWAREQ/S2qcF9XjUWrkShg6FN990k6smTYJbb4WkJN+VRR0FsIiIB1ndskhNSW10GDo1JZWsblkeqgpSXh6MHAkVFZCaCjNnwpln+q4qamkIWkTEg+SkZPKy80hNSd3h/u2zoGNqKVJNDYwf77aPrKiAs89220kqfH+UOmAREU8yOmVQOr6UQHGAkvUlsbkOeOlSd25vURG0bAl33eWWG+2l/m53FMAiIh4lJyWTc0yO7zKazlqYMcMdGVhZCWlpbpZzRhytXw4z/YgiIiJNs2mTOzpw5EgXvjk5sHChwreJ1AGLiMiey8+H7Gz4/HO3nvfRR2H48ITdTjIY6oBFRGT36uvh3nvhhBNc+B57LCxYACNGKHybSQEsIiI/bs0aOO88uPFGqKuDa65xZ/h27eq7spi22wA2xvzVGLPGGPPx9+7bzxjzmjHms4bf24e3TBER8eI//4GePeHf/4YOHeDFF+HBB91pRhKUPemAnwbO+cF9k4DXrbVHAK833BYRkXhRW+t2sTrrLFi92m0rWVQEmZm+K4sbuw1ga+3bwLof3D0AmNHw9QzgghDXJSIivpSWwkknwd13u+u7t94Kr78OnTv7riyuNHcWdKq19isAa+1XxpgDdvVEY8xYYCzAoYce2sy3ExGRiPjb32DsWLfU6JBDYM4c6NvXd1VxKeyTsKy106y16dba9I4dO4b77UREpDkqK91ZvYMHu/DNynJDzgrfsGluAJcbYw4CaPh9TehKEhGRiProI7eJxpNPuslVjz0Gzz0H++3nu7K41twAfhEY0fD1COCF0JQjIiIRY63bSOO442DJEujeHebPhyuu0NreCNiTZUi5wDygqzGmzBhzKXAXcKYx5jPgzIbbIiISK9atg1//2u3lvGULjBnjdrnq0cN3ZQljt5OwrLXZu3jo9BDXIiIikfDOOzBkCKxYAe3awRNPwMUX+64q4WgnLBGRRLFtG/z+93DqqS58+/RxE60Uvl7oMAYRkURQVgZDh8Jbb7nru5MmufW9SUm+K0tYCmARkXiXlweXXOKu+x54IMycCWec4buqhKchaBGReFVT4w5OOP98F77nnguLFil8o4Q6YBGReLR0KQwa5AI3KQnuvBOuvRb2Ut8VLRTAIiLxxFp4+mm3vKiqCn72M5g7F9LTfVcmP6AfhURE4sWmTW550ahRLnyHDIHCQoVvlFIHLCISD+bPh+xsKCmBlBR45BEYMUI7WkUxdcAiIrGsvh7uuQdOPNGFb69esGCBm/Ws8I1qCmARkVhVXg79+sHEiVBXB+PHw7x50LWr78pkD2gIWkQkFr36Kgwf7kK4Qwc38ap/f99VSROoAxYRiSW1tXDTTXD22S58Tz3VLTVS+MYcdcAiIrGipMRNtJo/363n/f3vYfJkaNHCd2XSDApgEZFYMHcuXHaZW2p0yCEwZw707eu7KgmChqBFRKJZZSWMHu06302bICvLnWCk8I156oBFRKLVokUweDAUF0Pr1vCnP7kuWMuL4oI6YBGRaGOt20ijd28Xvkcd5a77Xn65wjeOKIBFRKJJRYUbZr76atiyBcaOhfx8OOYY35VJiGkIWkQkWrz9ttu/uawMfvITeOIJuOgi31VJmKgDFhHxbds2t6TotNNc+PbpAwsXKnzjnDpgERGfyspc1/v22+767s03w9Sp7gxfiWsKYBERX158EUaOhHXr4MADYdYsOP1031VJhGgIWkQk0mpq3CSrAQNc+J57rltypPBNKOqARUQiqbjYre1dtMgNM999tzvFaC/1Q4lGASwiEgnWwlNPuc63qgoOP9xtL/nLX/quTDzRj1wiIuG2cSPk5MCll7rwHToUCgsVvglOHbCISDjNn++GnEtLISUFHnvMneMrCU8dsIhIONTXwx//CCee6MK3Vy/X9Sp8pYE6YBGRUCsvd0H76qvu9vjxbrJVq1Z+65KoogAWEQmlV1914VteDh06wNNPQ//+vquSKKQhaBGRUNi6FSZOhLPPduF76qluqZHCV3ZBHbCISLBKSiA72024atHCbSU5ebL7WmQXFMAiIsGYOxcuuww2bYJDD4U5c9zEK5Hd0BC0iEhzVFa6db3Z2S58L7wQiooUvrLH1AGLiDRVUZFb27t0KbRuDQ88AGPHutOMRPaQOmARkT1lLTzyiDuvd+lSOOooyM93Q9AKX2kiBbCIyJ6oqICsLLeX85YtLnTz8+Hoo31XJjFKQ9AiIrvz9tswZAiUlcFPfgLTp8PAgb6rkhinDlhEZFfq6tySotNOc+F7wgluba/CV0JAHbCISGNWrHBd7zvvuOu7U6a4MG6pb5sSGvqXJCLyQy+8AKNGwbp1cNBBMHMmnH6676okzmgIWkRku5oauOoquOACF779+rkhZ4WvhIE6YBERgCVL3Nrejz6CpCR3etH48bCX+hQJDwWwiCQ2a+Gpp9zyoqoqOPxwyM2F9HTflUmc0492IpK4Nm6EnBy3pWRVFQwbBoWFCl+JCHXAIpKY5s93Q86lpZCSAn/+swtgkQgJqgM2xlxrjPmfMeZjY0yuMaZ1qAoTEQmL+nr44x/doQmlpdCrl+t6Fb4SYc0OYGNMJ+AaIN1aezTQAhgcqsJEREJu9Wo491y46Sa3yca118K8eXDkkb4rkwQU7BB0SyDZGFMLtAFWBV+SiEgYvPIKDB8Oa9bA/vvDjBlumZGIJ83ugK21K4F7geXAV8BGa+2rP3yeMWasMabAGFOwdu3a5lcqItIcW7fCxIlwzjkufE87za3tVfiKZ8EMQbcHBgCHAQcDKcaYoT98nrV2mrU23Vqb3rFjx+ZXKiLSVJ9/Dn37wj33QIsWcPvt8NprcPDBvisTCWoS1hlAqbV2rbW2FngeOCE0ZYmIBCk3102wys+HQw91JxpNmeKCWCQKBBPAy4E+xpg2xhgDnA4sCU1ZIiLNVFnp9nHOyYFvvoELL4SiIneSkUgUafYkLGvth8aYfwCFQB2wEJgWqsJERJps0SIYNAiWLoXWreGBB2DsWHeakUiUCWoWtLX2d8DvQlSLiEjzWAuPPAI33OAmXf385zB3Lhx9tO/KRHZJW1GKSGyrqHCnF11zjQvfyy5zu1wpfCXKaStKEYldb70FQ4bAypWw774wfbq75isSA9QBi0jsqauD3/0OfvUrF74nnOAmWil8JYaoAxaR2LJihZvh/O67bnLVlCkwdSq01LcziS36FysiseP//s8tMVq/Hg46CGbNcl2wSAzSELSIRL/qahg3DrKyXPj26+eWHCl8JYapAxaR6LZkiVvbu3gxJCXB3XfDhAla2ysxTwEsItHJWnjySbe8qLoaDj/cre395S99VyYSEhqCFpHos3EjDB4MY8a48B02DAoLFb4SV9QBi0h0+fBDF75ffAH77AOPPeYCWCTOqAMWkehQX++u7/bt68L3F79wXa/CV+KUAlhE/Fu9Gs45ByZNcptsXHcdvP8+HHGE78pEwkZD0CLi1yuvwPDhsGYN7L8/zJjhlhmJxDl1wCLix9atcOONrvNds8at6V20SOErCUMdsIhE3uefQ3Y25OdDixZw220wcaL7WiRBKIBFJLJmz4YrroBvvoGf/hTmzHGHKYgkGA1Bi0hkbN4MI0fC0KEufAcOdCcYKXwlQakDFpHwKypy20l++im0bg0PPug22dB2kpLA1AGLSPhYCw89BL17u/D9+c+hoADGjlX4SsJTAItIeHz9NQwYAOPHuxnPl1/uJl39/Oe+KxOJChqCFpHQe/NNGDIEVq2CffeF6dPhwgt9VyUSVdQBi0jo1NXBb3/r1vSuWuUmWBUVKXxFGqEOWERCY/lyyMmB995z13enTIGpU6Glvs2INEb/M0QkeIEAjBoFGzbAwQfDrFlw2mm+qxKJahqCFpHmq66GK6+EX//ahW///m47SYWvyG6pAxaR5vnkE3du7+LFsPfe7ijB8eO1vEhkDymARaRprIUnn4RrrnEd8BFHwNy57vzeXaiqrSKwJEDphlLS2qeR1S2L5KTkCBYtEn0UwCKy5zZsgMsug7//3d0eMQIefhjatt3lH8lfmU9mbiblleXf3peakkpedh4ZnTLCXbFI1NI1YBHZMx98AL16ufDdZx+YOROefvpHw7e6tnqn8AUorywnMzeT6trqMBctEr0UwCLy4+rr4a67oG9f+OIL+OUvobDQHaqwG4HiwE7hu115ZTmB4kCIixWJHQpgEdm1r76Cs8+GyZNh2za47jq3zveII/boj5esLwnqcZF4pmvAItK4f/8bhg+HtWth//1hxgzo169JL5HWPi2ox0XimTpgEdnR1q1www1w7rkufH/1K7e2t4nhC5DVLYvUlNRGH0tNSSWrW1aw1YrELAWwiHxn2TI48US47z5o0QLuuANefdXtbtUMyUnJ5GXn7RTC22dBaymSJDINQYuIM3u2OzJw82b46U8hNxeOPz7ol83olEHp+FICxQFK1pf86DpgrReWRKIAFkl0mzfDVVe5a7wAAwfCE0+4YwRDJDkpmZxjcn70OVovLIlGQ9AiiWzhQresaMYMSE6GadPcOt8Qhu+e0HphSUQKYEkYVbVVzP5oNre/fTtzFs9J7G/q1sKDD0KfPvDpp3D00ZCfD2PGeNnLWeuFJRFpCFoSgoY3v+frr2HkSHjpJXf7iivcpKtkf9datV5YEpE6YIl7Gt78njffhJ49Xfjuuy889xw89pjX8AWtF5bEpACWuKfhTaCuDm65xa3pXbXKLTUqKnLn+EYBrReWRKQAlriX8MOby5fDqafC7be727fc4jrhn/7UZ1U70HphSUS6BixxL6GHN59/Hi691B0jePDBMGsWnHaa76oa1ZT1wiLxQAEscW/78GZjw9BxO7xZXQ3XXw9//rO73b8/PPWU29M5iu3JemGReKEhaIl7CTe8+ckn0Lu3C9+993bLjV58MerDVyTRqAOWhJAQw5vWwvTpMH6864CPPBLmzoVevXxXJiKNCCqAjTH7AtOBowELjLLWzgtFYSKhFtfDmxs2wNix8Oyz7vYll8DDD8M++3gtS0R2LdgO+EHg39bagcaYvYE2IahJRJpi3jzIzoYvv4S2bd3Q85AhvqsSkd1o9jVgY0w74GTgSQBr7VZr7YZQFSYiu7Ftmzsu8KSTXPimp0NhocJXJEYEMwkrDVgLPGWMWWiMmW6MSfnhk4wxY40xBcaYgrVr1wbxdiLyra++grPPhilTXBBffz289x4cfrjvykRkDwUTwC2BXwB/ttb2AiqBST98krV2mrU23Vqb3rFjxyDeTkQAePllt53k669Dx47wr3/Bvfe6Gc8iEjOCCeAyoMxa+2HD7X/gAllEwmHrVtfp9usHa9fCGWfAokVw7rm+KxORZmh2AFtrVwMrjDFdG+46HfgkJFWJyI6WLXP7N99/P7RoAXfeCa+8Agcd5LsyEWmmYGdBXw3MbpgBXQKMDL4kEdnBrFnuyMDNm6FLF8jNdef4ikhMCyqArbVFQHqIahGR79u8GcaNg2eecbcvugimTXPHCIpIzNNOWCLRqLAQBg+Gzz5zZ/U+9JA7VMEY35WJSIhoL2iRaGKt27v5+ONd+B59NBQUwOjRCl+ROKMAFokWX38N558PEya4Gc9XXAHz58NRR/muTETCQEPQItHgjTdg6FBYtcpd433ySfj1r31XJSJhpA5YxKe6OrjlFjj9dBe+J57o1vYqfEXinjpgEV++/BJycuD999313Vtugd/+Flrqv6VIItD/dBEfnnvOTazasAEOPhhmz4ZTT/VdlYhEkIagRSKpuhouvxwGDnTh27+/G3JW+IokHHXAIpHyv/+5tb0ff+wOTrjnHrj6ai0vEklQCmCRcLMWnnjCLS+qroYjj4S5c6FXL9+ViYhHGoIWCacNG2DQILjsMhe+l1wCCxYofEVEHbBI2MybB9nZbrbzPvvA44/DkCG+q/pWVW0VgSUBSjeUktY+jaxuWSQnJfsuSyRhKIBFQm3bNrj7brekaNs2SE93Jxgdfrjvyr6VvzKfzNxMyivLv70vNSWVvOw8MjpleKxMJHFoCFoklFatgrPOgilTXPjecAO8915UhW91bfVO4QtQXllOZm4m1bXVnioTSSwKYJFQ+de/oGdP+O9/oWNHePllN9N57719V7aDQHFgp/DdrryynEBxIMIViSQmDUGLBGvrVpg8Ge6/390+4wyYORMOPDBsbxnM9duS9SVBPS4ioaEAlrgS8YlFn33mJlotWOC2kLz9drjxRtgrfINLwV6/TWufFtTjIhIaCmCJGxGfWDRrljsycPNm6NLFTbTq0yf07/M9u7t+Wzq+dLc/cGR1yyI1JbXRYejUlFSyumWFtGYRaZyuAUtciOjEos2bYcQIGDbMfT1oEBQVhT18ITTXb5OTksnLziM1JXWH+7f/sKKlSCKRoQ5Y4sKeBFPOMTnBv1FhodtO8rPPIDkZHn4YRo2K2HaSobp+m9Epg9LxpQSKA5SsL9E6YBEPFMASF8I+schaeOghmDjRTbrq0cNtJ9m9e3Cv20ShvH6bnJQcmh9KRKRZNAQtcSGsE4vWroXMTLeX89atMG4cfPhhxMMXvrt+2xhdvxWJLQpgiQthC6Y33nBre//5T2jfHgIBeOQRaN06iGqbT9dvReKHhqAlLmwPpl3Ngm5yMNXVwdSpcMcdbvi5b1+YMwcOOSS0hTeDrt+KxAdjrY3Ym6Wnp9uCgoKIvZ8knura6uCD6csvIScH3n/free95Rb4zW/cOl8RkSYwxiyw1qY39pi+o0hcCXpi0XPPwejR7hjBTp1g9mw45ZTQFSgi0kDXgEXAndV7+eUwcKAL3/PPh0WLFL4iEjbqgEU+/tit7f3f/9zBCffd52Y6R2htr4gkJgWwJC5rYdo0t7yopga6dnVre4891ndlIpIANAQtiWn9erjoIjfsXFMDI0e6AxUUviISIeqAJfG8/747wWj5cmjbFv7yF3dbRCSC1AFL4ti2Df7wBzj5ZBe+GRmwcKHCV0S8UAcsiWHVKhg61O1sBXDDDS6M997bb10ikrAUwBL//vlPuOQS+PprOOAAeOYZOPts31WJSILTELTEry1b4NproX9/F75nnOHW9ip8RSQKqAOW+PTZZ25tb2Gh20Ly9tvhxhvd1pIiIlFAASzxZ+ZMuPJK2LwZunRxa3t79/ZdlYjIDtQOSPz45hsYPtz92rwZBg2CoiKFr4hEJXXAEh8KC13gLlsGycnw8MMwapS2kxSRqKUAll2qqq0isCRA6YbS6D1z1lp44AG46SaorYUePdyQc/fuvisTEflRCmBpVP7K/F0ebp/RKcNjZd+zdq1bXvSvf7nb48bBvfdC69ZeyxIR2RO6Biw7qa6t3il8Acory8nMzaS6ttpTZd/z3/9Cz54ufNu3h0AAHnlE4SsiMUMBLDsJFAd2Ct/tyivLCRQHIlzR99TWwpQpbk3vV1/BSSe5tb0XXOCvJhGRZtAQtOykZH1JUI+HzRdfQE4OzJvn1vPecov71VL/jEUk9ug7l+wkrX1aUI+HxT/+AaNHw8aN0KkTzJ4Np5wS+TpEREJEQ9Cyk6xuWaSmpDb6WGpKKlndsiJXTFUVXHaZO7t340bIzHRDzgpfEYlxQQewMaaFMWahMealUBQk/iUnJZOXnbdTCG+fBR2xpUgffwzHHQfTprlTix56CF54ATp0iMz7i4iEUSiGoMcDS4B2IXgtiRIZnTIoHV9KoDhAyfqSyK4DttaF7oQJUFMDXbu6tb3HHhv+9xYRiZCgAtgY0xk4D/gDcF1IKpKokZyUTM4xOZF90/XrYexYd80X3G5WDz0EKSmRrUNEJMyC7YAfACYCbXf1BGPMWGAswKGHHhrk20lce+89N8t5+XJo2xb+8hfIzvZdlYhIWDT7GrAxpj+wxlq74MeeZ62dZq1Nt9amd+zYsblvJ/Fs2zb4wx/cxKrly91136Iiha+IxLVgOuATgfONMf2A1kA7Y8wsa+3Q0JQmCWHVKhg6FN54w92eOBFuu81NuhIRiWPN7oCttZOttZ2ttV2AwcB/Fb7SJC+95A5PeOMNOOAAeOUVuPtuha+IJAStA5bI27LFzXDOzISKCjjzTLe296yzfFcmIhIxIdkJy1r7JvBmKF5L4tynn7pru4WFbgvJO+6A6693W0uKiCQQbUUpkfPMM3DllVBZCWlpkJvrJlyJiCQgtR0Sft98A8OGwYgRLnyzs2HhQoWviCQ0dcASXgsWwODBsGwZtGnjzuy95BIwxndlIiJeqQOW8Kivh/vvh+OPd+Hbo4cL45EjFb4iIiiAJRzWrIH+/d3kqtpauOoq+PBD6NbNd2UiIlFDQ9ASWq+/7jbWWL0a9tsP/vpXGDDAd1UiIlFHHbCERm0t3HyzW9O7ejWcfLJb26vwFRFplDpgCd4XX7hDFObNc+t5f/c7+M1voEUL35WJiEQtBbAE59lnYcwY2LgROneG2bNd9ysiIj9KQ9DSPFVV7tzeiy924TtggDvBSOErIrJH1AFL0y1e7Nb2fvIJtGoF994L48ZpeZGISBMogGXPWQuPPw7XXQc1NW5Z0dy50LOn78pERGKOhqBlz6xfDwMHur2ca2rg0kuhoEDhKyLSTOqAZffee8/Ncl6+HNq2hWnT3BC0iIg0mzpg2bVt2+D22+GUU1z4Hnecm2il8BURCZo6YGncypVuR6s333S3J050YZyU5LUsEZF4oQCWnb30kjuxqKICUlNh5ky3w5WIiISMhqDlO1u2wIQJkJnpwvess9x2kgpfEZGQUwcszqefumu7CxdCy5Zwxx3uNKO99DOaiEg4KIATnbXwzDNuI43KSjjsMMjNhd69fVcmIhLX1N4ksm++gWHD3PXeykrIznYdsMJXRCTs1AEnqgUL3JDzsmXQpg088ogLYm0nKSISEQrgRFNfD3/6E0ye7M7w7dkT5s6l6meHElg8h9INpaS1TyOrWxbJScm+qxURiVsK4ESyZo3rcl9+2d2++mr44x/Jr1hM5oNplFeWf/vU1JRU8sfYbK0AAA3FSURBVLLzyOiU4adWEZE4p2vAieL11123+/LLsN9+8MIL8NBDVLewZOZm7hC+AOWV5WTmZlJdW+2pYBGR+KYAjne1tXDzzW4t7+rV7rzeRYvg/PMBCBQHdgrf7corywkUByJZrYhIwtAQdDz74gt3iMK8eW4979SpMGUKtGjx7VNK1pf86Evs7nEREWkeBXC8evZZGDMGNm6Ezp1hzhw46aSdnpbWPu1HX2Z3j4uISPNoCDreVFW54L34Yhe+F1zghpwbCV+ArG5ZpKakNvpYakoqWd2ywlmtiEjCUgDHk8WLIT0dpk+HVq3g0Ufh+efdpKtdSE5KJi87b6cQ3j4LWkuRRETCQ0PQ8cBaePxxuPZad6BCt27wt79Bjx579MczOmVQOr6UQHGAkvUlWgcsIhIBCuBYt24djB4NgYbZyqNHwwMPQEpKk14mOSmZnGNywlCgiIg0RgEcy959181yXrEC2rWDadNg0CDfVYmIyB7QNeBYtG0b3HYbnHKKC9/evaGoSOErIhJD1AHHmpUrYehQePNNd3DCpElw662QlOS7MhERaQIFcCzJy4ORI6GiAlJTYeZMt8OViIjEHA1Bx4ItW2D8eLd9ZEUFnH22W9ur8BURiVnqgKPd0qXu3N6iImjZEu68E667zm0tKSIiMUsBHK2shWeegXHjoLIS0tJg7lzI0PGAIiLxQG1UNNq0yU20uuQSF77Z2bBwocJXRCSOqAOONvn5LnA//xzatHHbSY4Y4WY8i4hI3FAHHC3q6+G+++CEE1z4HnssFBa6LljhKyISdxTA0aC8HM47D264Aerq4Oqr3Rm+Xbv6rkxERMJEQ9C+vfYaDBvmQrhDB3jqKcjM9F2ViIiEmTpgX2prYfJkt6a3vNxtK7lokcJXRCRBNLsDNsYcAjwDHAjUA9OstQ+GqrC4VlrqDlH44AO3nvf3v4ebb4YWLbyVVFVbRWBJgNINpTqOUEQkAoIZgq4DrrfWFhpj2gILjDGvWWs/CVFt8envf4cxY9xSo0MOgTlzoG9fryXlr8wnMzeT8sryb+9LTUklLzuPjE5a+iQiEg7NHoK21n5lrS1s+PobYAnQKVSFxZ2qKhe8gwa58M3KcrtbeQ7f6trqncIXoLyynMzcTKprqz1VJiIS30JyDdgY0wXoBXwYiteLOx99BOnpMH06tGoFjz0Gzz0H++3nuzICxYGdwne78spyAsWBCFckIpIYgg5gY8w+wHPABGvtpkYeH2uMKTDGFKxduzbYt4st1rqwPe44WLIEuneH+fPhiiuiZm1vyfqSoB4XEZHmCSqAjTFJuPCdba19vrHnWGunWWvTrbXpHTt2DObtYsu6dXDhhW4v5y1bYPRot8tVjx6+K9tBWvu0oB4XEZHmaXYAG2MM8CSwxFp7f+hKigPvvON2sgoEoF07d4jCE09ASorvynaS1S2L1JTURh9LTUklq1tWhCsSEUkMwXTAJwLDgF8ZY4oafvULUV2xads2uPVWOPVUWLECevd2E60GDfJd2S4lJyWTl523UwhvnwWtpUgiIuHR7GVI1tp3gei4kBkNysrcCUZvveWu706a5MI4Kcl3ZbuV0SmD0vGlBIoDlKwv0TpgEZEI0FaUofDiizBypLvum5oKM2fCmWf6rqpJkpOSyTkmx3cZIiIJQ1tRBqOmBq65BgYMcOF7zjluyVGMha+IiESeOuDmWroUBg9213iTkuDOO+Haa93WkiIiIruhAG4qa2HGDLjqKqishJ/9zM1yTk/3XZmIiMQQBXBTbNrkNtGYM8fdHjLEbbTRrp3fukREolBtbS1lZWXU1NT4LiXsWrduTefOnUlqwsRbBfCeys+H7Gz4/HO3nvfRR2H48KjZ0UpEJNqUlZXRtm1bunTpgonj75XWWioqKigrK+Owww7b4z+nC5a7U18P994LJ5zgwvfYY2HBAhgxQuErIvIjampq6NChQ1yHL4Axhg4dOjS501cA/5jycujXD268EerqYPx4d4Zv166+KxMRiQnxHr7bNefvqSHoXXntNRg2zIVwhw7w9NPQv7/vqkREJE6oA/6h2lq3i9VZZ7nwPeUUWLRI4SsiIiGlDvj7SkvdRKsPP3TreadOhZtvhhYtfFcmIiLNMHXqVD744ANatnRxV1dXR58+fRq9D2jS/VOnTg2qNgXwdn/7G4wd65YaHXKIW2rUt6/vqkRE4kO4rgVbu9unzJ07l3333ReADRs28MADDzR6366e+2P3B0ND0JWV7qzewYNd+F5wgdvdSuErIiJhlNgd8EcfuaMCi4uhVSv405/g8su1vEhEJNT2oFNNNInZAVvrNtI47jgXvt27u402rrhC4SsiIhGReAG8bh1kZbm9nLdsgTFjoKAAjjnGd2UiIpJAEmsI+p13ICcHysrgJz+BJ56Aiy7yXZWIiCSgxOiAt22DW2+FU0914dunj5topfAVERFP4r8DLiuDoUPhrbfc9d1Jk1wYN+HEChERiU0HHHAAw4cPZ6+Gs9rr6+s555xzGr0PaPL9wTA2gjPT0tPTbUFBQcTejxdfhJEj3XXfAw+EmTPhjDMi9/4iIglsyZIldO/e3XcZEdPY39cYs8Ba2+iB8fE5BF1TA9dcAwMGuPA991y3naTCV0REokT8DUEXF7tNNRYtcsPMd90FEya4rSVFRESiRPwEsLXw1FNw9dVQVQWHHw65uZDeaOcvIiLiVXy0hRs3wpAhcOmlLnyHDoXCQoWviIhErdjvgOfPd0POpaWQkuJ2uBoxwndVIiLSRFW1VQSWBCjdUEpa+zSyumWRnJTsu6ywid0Arq+H++5zxwXW1UGvXjB3Lhx5pO/KRESkifJX5pOZm0l5Zfm396WmpJKXnUdGp4xmv66OIwy1LVvcDOdXXnG3J0xwk61atfJbl4iINFl1bfVO4QtQXllOZm4mpeNLg+qEdRxhKLVqBQcfDB06QF6eO8VI4SsiEpMCxYGdwne78spyAsWBCFcUGbEZwAAPP+yWGvXv77sSEREJQsn6kqAej1WxOQQNbsJVSorvKkREJEhp7dOCejxWxW4HLCIicSGrWxapKamNPpaakkpWt6wIVxQZCmAREfEqOSmZvOy8nUJ4+yzoeF2KFLtD0CIiEjcyOmVQOr6UQHGAkvUlWgcsIiISKclJyeQckxPS19RxhA0ifhyhiIh4o+MIE/E4QhERiQqRbPJ8as7fUwEsIiJh0bp1ayoqKuI+hK21VFRU0Lp16yb9OV0DFhGRsOjcuTNlZWWsXbvWdylh17p1azp37tykP6MAFhGRsEhKSuKwww7zXUbU0hC0iIiIBwpgERERDxTAIiIiHkR0HbAxZi3wZQhfcn/g6xC+XqzT57EjfR7f0WexI30e39FnsaNQfx4/tdZ2bOyBiAZwqBljCna1wDkR6fPYkT6P7+iz2JE+j+/os9hRJD8PDUGLiIh4oAAWERHxINYDeJrvAqKMPo8d6fP4jj6LHenz+I4+ix1F7POI6WvAIiIisSrWO2AREZGYFLMBbIw5xxiz1BizzBgzyXc9vhhjDjHGvGGMWWKM+Z8xZrzvmqKBMaaFMWahMeYl37X4ZozZ1xjzD2NMccO/k+N91+SLMebahv8nHxtjco0xTds9P8YZY/5qjFljjPn4e/ftZ4x5zRjzWcPv7X3WGEm7+Dzuafi/8pExJmCM2Tdc7x+TAWyMaQE8CpwLHAVkG2OO8luVN3XA9dba7kAfYFwCfxbfNx5Y4ruIKPEg8G9rbTegJwn6uRhjOgHXAOnW2qOBFsBgv1VF3NPAD0+SnwS8bq09Ani94XaieJqdP4/XgKOttT2AT4HJ4XrzmAxg4DhgmbW2xFq7FZgLDPBckxfW2q+stYUNX3+D++bayW9VfhljOgPnAdN91+KbMaYdcDLwJIC1dqu1doPfqrxqCSQbY1oCbYBVnuuJKGvt28C6H9w9AJjR8PUM4IKIFuVRY5+HtfZVa21dw80PgKYdcdQEsRrAnYAV37tdRoKHDoAxpgvQC/jQbyXePQBMBOp9FxIF0oC1wFMNQ/LTjTEpvovywVq7ErgXWA58BWy01r7qt6qokGqt/QrcD/TAAZ7riSajgJfD9eKxGsCmkfsSejq3MWYf4DlggrV2k+96fDHG9AfWWGsX+K4lSrQEfgH82VrbC6gksYYYv9VwbXMAcBhwMJBijBnqtyqJVsaYKbhLfLPD9R6xGsBlwCHfu92ZBBtK+j5jTBIufGdba5/3XY9nJwLnG2O+wF2a+JUxZpbfkrwqA8qstdtHRf6BC+REdAZQaq1da62tBZ4HTvBcUzQoN8YcBNDw+xrP9XhnjBkB9AeG2DCu1Y3VAM4HjjDGHGaM2Rs3keJFzzV5YYwxuOt7S6y19/uuxzdr7WRrbWdrbRfcv4v/WmsTtsux1q4GVhhjujbcdTrwiceSfFoO9DHGtGn4f3M6CToh7QdeBEY0fD0CeMFjLd4ZY84BbgLOt9ZWhfO9YjKAGy6QXwW8gvsP9Hdr7f/8VuXNicAwXKdX1PCrn++iJKpcDcw2xnwEHAvc4bkeLxpGAf4BFAKLcd//EmoXKGNMLjAP6GqMKTPGXArcBZxpjPkMOLPhdkLYxefxCNAWeK3h++njYXt/7YQlIiISeTHZAYuIiMQ6BbCIiIgHCmAREREPFMAiIiIeKIBFREQ8UACLiIh4oAAWERHxQAEsIiLiwf8DtZqFT56wRI8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 引入最小二乘法算法的简单\n",
    "\n",
    "import numpy as np   ##科学计算库 \n",
    "import scipy as sp   ##在numpy基础上实现的部分算法库\n",
    "import matplotlib.pyplot as plt  ##绘图库\n",
    "from scipy.optimize import leastsq  ##引入最小二乘法算法\n",
    "\n",
    "\n",
    "##样本数据(Xi,Yi)，需要转换成数组(列表)形式\n",
    "Xi=np.array([6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2])\n",
    "Yi=np.array([5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3])\n",
    "\n",
    "'''\n",
    "    设定拟合函数和偏差函数\n",
    "    函数的形状确定过程：\n",
    "    1.先画样本图像\n",
    "    2.根据样本图像大致形状确定函数形式(直线、抛物线、正弦余弦等)\n",
    "'''\n",
    "\n",
    "##需要拟合的函数func :指定函数的形状\n",
    "def func(p,x):\n",
    "    k,b=p\n",
    "    return k*x+b\n",
    "\n",
    "##偏差函数：x,y都是列表:这里的x,y更上面的Xi,Yi中是一一对应的\n",
    "def error(p,x,y):\n",
    "    return func(p,x)-y\n",
    "\n",
    "'''\n",
    "    主要部分：附带部分说明\n",
    "    1.leastsq函数的返回值tuple，第一个元素是求解结果，第二个是求解的代价值(个人理解)\n",
    "    2.官网的原话（第二个值）：Value of the cost function at the solution\n",
    "    3.实例：Para=>(array([ 0.61349535,  1.79409255]), 3)\n",
    "    4.返回值元组中第一个值的数量跟需要求解的参数的数量一致\n",
    "'''\n",
    "\n",
    "#k,b的初始值，可以任意设定,经过几次试验，发现p0的值会影响cost的值：Para[1]\n",
    "p0=[1,20]\n",
    "\n",
    "#把error函数中除了p0以外的参数打包到args中(使用要求)\n",
    "Para=leastsq(error,p0,args=(Xi,Yi))\n",
    "\n",
    "#读取结果\n",
    "k,b=Para[0]\n",
    "print(\"k=\",k,\"b=\",b)\n",
    "print(\"cost：\"+str(Para[1]))\n",
    "print(\"求解的拟合直线为:\")\n",
    "print(\"y=\"+str(round(k,2))+\"x+\"+str(round(b,2)))\n",
    "\n",
    "#画样本点\n",
    "plt.figure(figsize=(8,6)) ##指定图像比例： 8：6\n",
    "plt.scatter(Xi,Yi,color=\"green\",label=\"样本数据\",linewidth=2) \n",
    "\n",
    "#画拟合直线\n",
    "x=np.linspace(0,12,100) ##在0-15直接画100个连续点\n",
    "y=k*x+b ##函数式\n",
    "plt.plot(x,y,color=\"red\",label=\"拟合直线\",linewidth=2) \n",
    "plt.legend(loc='lower right') #绘制图例\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
